// 文件名：mdspan_matrix_multiply.cpp
// 目的：使用 C++23 std::mdspan 为现有连续内存提供 2D 视图，并实现矩阵相乘 C = A * B。
// 核心概念：
//   - std::mdspan<Element, Extents[, Layout[, Accessor]]> 是对已有内存的“非拥有视图”；
//   - 本例使用动态维度（std::dynamic_extent）构建二维视图：A(m×k)、B(k×n)、C(m×n)；
//   - 缺省布局 layout_right（行主序），使线性内存以 row-major 方式索引。
// API/实现：
//   - matrix_view：只读视图，matrix_view_mutable：可写视图；
//   - matmul(A,B,C)：检查维度一致性，然后三重循环计算 C(i,j)；
// 注意事项：
//   - mdspan 不管理内存所有权，调用者需保证底层 storage 的生命周期；
//   - 若需要其他布局（列主序、跨距不规则），可通过 layout_left / layout_stride 适配；
//   - 性能优化（可选）：可考虑块状分块（tiling）、SIMD、并行等。
// 复杂度：O(M*K*N)。
#include <algorithm>
#include <iostream>
#include <mdspan>
#include <vector>

using matrix_view =
    std::mdspan<const double,
                std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>>;

using matrix_view_mutable =
    std::mdspan<double,
                std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>>;

void matmul(matrix_view A, matrix_view B, matrix_view_mutable C) {
    auto [m, k] = std::array{A.extent(0), A.extent(1)};
    auto [k2, n] = std::array{B.extent(0), B.extent(1)};
    if (k != k2 || C.extent(0) != m || C.extent(1) != n) {
        throw std::runtime_error("dimension mismatch");
    }

    for (std::size_t i = 0; i < m; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            double acc = 0.0;
            for (std::size_t t = 0; t < k; ++t) {
                acc += A(i, t) * B(t, j);
            }
            C(i, j) = acc;
        }
    }
}

int main() {
    constexpr std::size_t M = 2, K = 3, N = 2;
    std::vector<double> storageA{1, 2, 3,
                                 4, 5, 6}; // 行主序 (2x3): A(0,:) = [1 2 3], A(1,:) = [4 5 6]
    std::vector<double> storageB{7,  8,
                                 9,  10,
                                 11, 12}; // (3x2)
    std::vector<double> storageC(M * N, 0.0);

    matrix_view A(storageA.data(), M, K); // 默认 layout_right
    matrix_view B(storageB.data(), K, N);
    matrix_view_mutable C(storageC.data(), M, N);

    matmul(A, B, C);

    for (std::size_t i = 0; i < M; ++i) {
        for (std::size_t j = 0; j < N; ++j) {
            std::cout << C(i, j) << ' ';
        }
        std::cout << '\n';
    }
    return 0;
}